home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zuni1.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  10.0 KB  |  243 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((zeror 0.0) (zeroi 0.0) (coner 1.0))
  12.   (declare (type double-float coner zeroi zeror))
  13.   (defun zuni1 (zr zi fnu kode n yr yi nz nlast fnul tol elim alim)
  14.     (declare (type (simple-array double-float (*)) yr yi)
  15.              (type f2cl-lib:integer4 kode n nz nlast)
  16.              (type double-float zr zi fnu fnul tol elim alim))
  17.     (prog ((bry (make-array 3 :element-type 'double-float))
  18.            (cwrkr (make-array 16 :element-type 'double-float))
  19.            (cwrki (make-array 16 :element-type 'double-float))
  20.            (cssr (make-array 3 :element-type 'double-float))
  21.            (csrr (make-array 3 :element-type 'double-float))
  22.            (cyr (make-array 2 :element-type 'double-float))
  23.            (cyi (make-array 2 :element-type 'double-float)) (i 0) (iflag 0)
  24.            (init 0) (k 0) (m 0) (nd 0) (nn 0) (nuf 0) (nw 0) (aphi 0.0)
  25.            (ascle 0.0) (crsc 0.0) (cscl 0.0) (c1r 0.0) (c2i 0.0) (c2m 0.0)
  26.            (c2r 0.0) (fn 0.0) (phii 0.0) (phir 0.0) (rast 0.0) (rs1 0.0)
  27.            (rzi 0.0) (rzr 0.0) (sti 0.0) (str 0.0) (sumi 0.0) (sumr 0.0)
  28.            (s1i 0.0) (s1r 0.0) (s2i 0.0) (s2r 0.0) (zeta1i 0.0) (zeta1r 0.0)
  29.            (zeta2i 0.0) (zeta2r 0.0))
  30.       (declare (type (simple-array double-float (2)) cyi cyr)
  31.                (type (simple-array double-float (16)) cwrkr cwrki)
  32.                (type (simple-array double-float (3)) cssr csrr bry)
  33.                (type double-float zeta2r zeta2i zeta1r zeta1i s2r s2i s1r s1i
  34.                 sumr sumi str sti rzr rzi rs1 rast phir phii fn c2r c2m c2i c1r
  35.                 cscl crsc ascle aphi)
  36.                (type f2cl-lib:integer4 nw nuf nn nd m k init iflag i))
  37.       (setf nz 0)
  38.       (setf nd n)
  39.       (setf nlast 0)
  40.       (setf cscl (/ 1.0 tol))
  41.       (setf crsc tol)
  42.       (f2cl-lib:fset (f2cl-lib:fref cssr (1) ((1 3))) cscl)
  43.       (f2cl-lib:fset (f2cl-lib:fref cssr (2) ((1 3))) coner)
  44.       (f2cl-lib:fset (f2cl-lib:fref cssr (3) ((1 3))) crsc)
  45.       (f2cl-lib:fset (f2cl-lib:fref csrr (1) ((1 3))) crsc)
  46.       (f2cl-lib:fset (f2cl-lib:fref csrr (2) ((1 3))) coner)
  47.       (f2cl-lib:fset (f2cl-lib:fref csrr (3) ((1 3))) cscl)
  48.       (f2cl-lib:fset (f2cl-lib:fref bry (1) ((1 3)))
  49.                      (/ (* 1000.0 (f2cl-lib:d1mach 1)) tol))
  50.       (setf fn (max fnu 1.0))
  51.       (setf init 0)
  52.       (multiple-value-bind
  53.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  54.            var-11 var-12 var-13 var-14 var-15 var-16)
  55.           (zunik zr zi fn 1 1 tol init phir phii zeta1r zeta1i zeta2r zeta2i
  56.            sumr sumi cwrkr cwrki)
  57.         (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-15 var-16))
  58.         (setf init var-6)
  59.         (setf phir var-7)
  60.         (setf phii var-8)
  61.         (setf zeta1r var-9)
  62.         (setf zeta1i var-10)
  63.         (setf zeta2r var-11)
  64.         (setf zeta2i var-12)
  65.         (setf sumr var-13)
  66.         (setf sumi var-14))
  67.       (if (= kode 1) (go label10))
  68.       (setf str (+ zr zeta2r))
  69.       (setf sti (+ zi zeta2i))
  70.       (setf rast (/ fn (zabs str sti)))
  71.       (setf str (* str rast rast))
  72.       (setf sti (* (- sti) rast rast))
  73.       (setf s1r (- str zeta1r))
  74.       (setf s1i (- sti zeta1i))
  75.       (go label20)
  76.      label10
  77.       (setf s1r (- zeta2r zeta1r))
  78.       (setf s1i (- zeta2i zeta1i))
  79.      label20
  80.       (setf rs1 s1r)
  81.       (if (> (abs rs1) elim) (go label130))
  82.      label30
  83.       (setf nn (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 nd)))
  84.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  85.                     ((> i nn) nil)
  86.         (tagbody
  87.           (setf fn (+ fnu (f2cl-lib:int-sub nd i)))
  88.           (setf init 0)
  89.           (multiple-value-bind
  90.               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
  91.                var-10 var-11 var-12 var-13 var-14 var-15 var-16)
  92.               (zunik zr zi fn 1 0 tol init phir phii zeta1r zeta1i zeta2r
  93.                zeta2i sumr sumi cwrkr cwrki)
  94.             (declare
  95.              (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-15 var-16))
  96.             (setf init var-6)
  97.             (setf phir var-7)
  98.             (setf phii var-8)
  99.             (setf zeta1r var-9)
  100.             (setf zeta1i var-10)
  101.             (setf zeta2r var-11)
  102.             (setf zeta2i var-12)
  103.             (setf sumr var-13)
  104.             (setf sumi var-14))
  105.           (if (= kode 1) (go label40))
  106.           (setf str (+ zr zeta2r))
  107.           (setf sti (+ zi zeta2i))
  108.           (setf rast (/ fn (zabs str sti)))
  109.           (setf str (* str rast rast))
  110.           (setf sti (* (- sti) rast rast))
  111.           (setf s1r (- str zeta1r))
  112.           (setf s1i (+ (- sti zeta1i) zi))
  113.           (go label50)
  114.          label40
  115.           (setf s1r (- zeta2r zeta1r))
  116.           (setf s1i (- zeta2i zeta1i))
  117.          label50
  118.           (setf rs1 s1r)
  119.           (if (> (abs rs1) elim) (go label110))
  120.           (if (= i 1) (setf iflag 2))
  121.           (if (< (abs rs1) alim) (go label60))
  122.           (setf aphi (zabs phir phii))
  123.           (setf rs1 (+ rs1 (f2cl-lib:flog aphi)))
  124.           (if (> (abs rs1) elim) (go label110))
  125.           (if (= i 1) (setf iflag 1))
  126.           (if (< rs1 0.0) (go label60))
  127.           (if (= i 1) (setf iflag 3))
  128.          label60
  129.           (setf s2r (- (* phir sumr) (* phii sumi)))
  130.           (setf s2i (+ (* phir sumi) (* phii sumr)))
  131.           (setf str (* (exp s1r) (f2cl-lib:fref cssr (iflag) ((1 3)))))
  132.           (setf s1r (* str (cos s1i)))
  133.           (setf s1i (* str (sin s1i)))
  134.           (setf str (- (* s2r s1r) (* s2i s1i)))
  135.           (setf s2i (+ (* s2r s1i) (* s2i s1r)))
  136.           (setf s2r str)
  137.           (if (/= iflag 1) (go label70))
  138.           (multiple-value-bind
  139.               (var-0 var-1 var-2 var-3 var-4)
  140.               (zuchk s2r s2i nw (f2cl-lib:fref bry (1) ((1 3))) tol)
  141.             (declare (ignore var-0 var-1 var-3 var-4))
  142.             (setf nw var-2))
  143.           (if (/= nw 0) (go label110))
  144.          label70
  145.           (f2cl-lib:fset (f2cl-lib:fref cyr (i) ((1 2))) s2r)
  146.           (f2cl-lib:fset (f2cl-lib:fref cyi (i) ((1 2))) s2i)
  147.           (setf m (f2cl-lib:int-add (f2cl-lib:int-sub nd i) 1))
  148.           (f2cl-lib:fset (f2cl-lib:fref yr (m) ((1 n)))
  149.                          (* s2r (f2cl-lib:fref csrr (iflag) ((1 3)))))
  150.           (f2cl-lib:fset (f2cl-lib:fref yi (m) ((1 n)))
  151.                          (* s2i (f2cl-lib:fref csrr (iflag) ((1 3)))))
  152.          label80))
  153.       (if (<= nd 2) (go label100))
  154.       (setf rast (/ 1.0 (zabs zr zi)))
  155.       (setf str (* zr rast))
  156.       (setf sti (* (- zi) rast))
  157.       (setf rzr (* (+ str str) rast))
  158.       (setf rzi (* (+ sti sti) rast))
  159.       (f2cl-lib:fset (f2cl-lib:fref bry (2) ((1 3)))
  160.                      (/ 1.0 (f2cl-lib:fref bry (1) ((1 3)))))
  161.       (f2cl-lib:fset (f2cl-lib:fref bry (3) ((1 3))) (f2cl-lib:d1mach 2))
  162.       (setf s1r (f2cl-lib:fref cyr (1) ((1 2))))
  163.       (setf s1i (f2cl-lib:fref cyi (1) ((1 2))))
  164.       (setf s2r (f2cl-lib:fref cyr (2) ((1 2))))
  165.       (setf s2i (f2cl-lib:fref cyi (2) ((1 2))))
  166.       (setf c1r (f2cl-lib:fref csrr (iflag) ((1 3))))
  167.       (setf ascle (f2cl-lib:fref bry (iflag) ((1 3))))
  168.       (setf k (f2cl-lib:int-sub nd 2))
  169.       (setf fn (coerce (the f2cl-lib:integer4 k) 'double-float))
  170.       (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
  171.                     ((> i nd) nil)
  172.         (tagbody
  173.           (setf c2r s2r)
  174.           (setf c2i s2i)
  175.           (setf s2r (+ s1r (* (+ fnu fn) (- (* rzr c2r) (* rzi c2i)))))
  176.           (setf s2i (+ s1i (* (+ fnu fn) (+ (* rzr c2i) (* rzi c2r)))))
  177.           (setf s1r c2r)
  178.           (setf s1i c2i)
  179.           (setf c2r (* s2r c1r))
  180.           (setf c2i (* s2i c1r))
  181.           (f2cl-lib:fset (f2cl-lib:fref yr (k) ((1 n))) c2r)
  182.           (f2cl-lib:fset (f2cl-lib:fref yi (k) ((1 n))) c2i)
  183.           (setf k (f2cl-lib:int-sub k 1))
  184.           (setf fn (- fn 1.0))
  185.           (if (>= iflag 3) (go label90))
  186.           (setf str (coerce (abs c2r) 'double-float))
  187.           (setf sti (coerce (abs c2i) 'double-float))
  188.           (setf c2m (max str sti))
  189.           (if (<= c2m ascle) (go label90))
  190.           (setf iflag (f2cl-lib:int-add iflag 1))
  191.           (setf ascle (f2cl-lib:fref bry (iflag) ((1 3))))
  192.           (setf s1r (* s1r c1r))
  193.           (setf s1i (* s1i c1r))
  194.           (setf s2r c2r)
  195.           (setf s2i c2i)
  196.           (setf s1r (* s1r (f2cl-lib:fref cssr (iflag) ((1 3)))))
  197.           (setf s1i (* s1i (f2cl-lib:fref cssr (iflag) ((1 3)))))
  198.           (setf s2r (* s2r (f2cl-lib:fref cssr (iflag) ((1 3)))))
  199.           (setf s2i (* s2i (f2cl-lib:fref cssr (iflag) ((1 3)))))
  200.           (setf c1r (f2cl-lib:fref csrr (iflag) ((1 3))))
  201.          label90))
  202.      label100
  203.       (go end_label)
  204.      label110
  205.       (if (> rs1 0.0) (go label120))
  206.       (f2cl-lib:fset (f2cl-lib:fref yr (nd) ((1 n))) zeror)
  207.       (f2cl-lib:fset (f2cl-lib:fref yi (nd) ((1 n))) zeroi)
  208.       (setf nz (f2cl-lib:int-add nz 1))
  209.       (setf nd (f2cl-lib:int-sub nd 1))
  210.       (if (= nd 0) (go label100))
  211.       (multiple-value-bind
  212.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  213.            var-11)
  214.           (zuoik zr zi fnu kode 1 nd yr yi nuf tol elim alim)
  215.         (declare
  216.          (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-9 var-10
  217.           var-11))
  218.         (setf nuf var-8))
  219.       (if (< nuf 0) (go label120))
  220.       (setf nd (f2cl-lib:int-sub nd nuf))
  221.       (setf nz (f2cl-lib:int-add nz nuf))
  222.       (if (= nd 0) (go label100))
  223.       (setf fn (+ fnu (f2cl-lib:int-sub nd 1)))
  224.       (if (>= fn fnul) (go label30))
  225.       (setf nlast nd)
  226.       (go end_label)
  227.      label120
  228.       (setf nz -1)
  229.       (go end_label)
  230.      label130
  231.       (if (> rs1 0.0) (go label120))
  232.       (setf nz n)
  233.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  234.                     ((> i n) nil)
  235.         (tagbody
  236.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n))) zeror)
  237.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n))) zeroi)
  238.          label140))
  239.       (go end_label)
  240.      end_label
  241.       (return (values nil nil nil nil nil nil nil nz nlast nil nil nil nil)))))
  242.  
  243.